#include "matrix3D.cuh"
#include "timer.cuh"
#include <cstdio>

__global__ void gaus(short* a, int cur_dim_size, int step, char dim_number)
{
    const unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
    int align;

    if (dim_number == 0)
        align = cur_dim_size * tid + 1;

    else if (dim_number == 1)
        align = step * cur_dim_size * (tid/step) + step + (tid % step);
    else
        align = tid + step;

    for (int i=1; i < cur_dim_size - 1; i++)
        a[align + step * (i - 1)] = a[align + step * (i - 2)] * 0.25 + a[align + step * (i - 1)] *0.5 + a[align + step * i] * 0.25;
}

void controlSmooth(short* input, short* output, int dimx, int dimy, int dimz)
{
    short* buffer;
    cudaMalloc((void**)&buffer, dimx * dimy * dimz * 2);
    cudaMemcpy(buffer,input,dimx*dimy*dimz*2,cudaMemcpyHostToDevice);

    Timer t;
    cudaEvent_t start = t.go();

    if(((dimy*dimz) % 64)==0)
        gaus<<<(dimy*dimz)/64,64>>>(buffer,dimx,1,0);
    else
        gaus<<<(int)((dimy*dimz)/64)+1,64>>>(buffer,dimx,1,0);



    if(((dimx*dimz) % 64)==0)
        gaus<<<(dimx*dimz)/64,64>>>(buffer,dimy,dimx,1);
    else
        gaus<<<(int)((dimx*dimz)/64)+1,64>>>(buffer,dimy,dimx,1);



    if(((dimx*dimz) % 64)==0)
        gaus<<<(dimx*dimz)/64,64>>>(buffer,dimz,dimx*dimy,2);
    else
        gaus<<<(int)((dimx*dimz)/64)+1,64>>>(buffer,dimz,dimx*dimy,2);

    printf("Gaussian filtration finished after : %f Ms \n", t.get(start));

    cudaMemcpy(output, buffer, dimx * dimy * dimz * 2,cudaMemcpyDeviceToHost);
    cudaFree(buffer);
}
